Para replicar las secciones [1.] OpenStreetMap y
[2.] Otras fuentes de datos, debe descargar el
siguiente proyecto
de R y abrir el archivo clase-03.Rproj.
OpenStreetMap (OSM) es un proyecto de mapeo de acceso abierto global, gratuito y licenciado bajo ODbL Licence. OSM recopila información geográfica capturada con dispositivos GPS móviles, ortofotografías y otras fuentes de información libre.
## Llamar pacman (contiene la función p_load)
require(pacman)
## Llama/instala-llama las librerías listadas
p_load(tidyverse,rio,skimr,
viridis, # paletas de colores
sf, # Leer/escribir/manipular datos espaciales
leaflet, # Visualizaciones dinámicas
tmaptools, # geocode_OSM()
osmdata) # Get OSM's dataLa función geocode_OSM() de la librería
tmaptools se conecta a la API de OpenStreetMap
y retorna el vértice con la coordenada geográfica del sitio/dirección
buscado.
## Buscar un lugar público por el nombre
geocode_OSM("Casa de Nariño, Bogotá")## $query
## [1] "Casa de Nariño, Bogotá"
##
## $coords
## x y
## -74.077549 4.595506
##
## $bbox
## xmin ymin xmax ymax
## -74.077962 4.595103 -74.077078 4.595872
Puede adicionar el argumento as.sf=T para generar un
objeto de clase sf:
## geocode_OSM no reconoce el caracter #, en su lugar se usa %23%
cbd <- geocode_OSM("Centro Internacional, Bogotá", as.sf=T)
cbd## Simple feature collection with 1 feature and 7 fields
## Active geometry column: point
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: -74.07057 ymin: 4.613615 xmax: -74.07057 ymax: 4.613615
## Geodetic CRS: WGS 84
## query lat lon lat_min lat_max lon_min
## 1 Centro Internacional, Bogotá 4.613615 -74.07057 4.613373 4.613889 -74.07098
## lon_max bbox point
## 1 -74.07023 POLYGON ((-74.07098 4.61337... POINT (-74.07057 4.613615)
Puede visualizar el objeto point usando la librería
leaflet:
## la función addTiles adiciona la capa de OpenStreetMap
leaflet() %>% addTiles() %>% addCircles(data=cbd)osmdataosmdata es una librería de R que permite descargar y
usar datos de OpenStreetMap (OSM).
Puede acceder a la lista de features disponibles en
OSM aquí. En R
puede obtener un vector con los nombres de los features
usando la función available_features():
available_features() %>% head(20)## [1] "4wd_only" "abandoned"
## [3] "abutters" "access"
## [5] "addr" "addr:city"
## [7] "addr:conscriptionnumber" "addr:country"
## [9] "addr:county" "addr:district"
## [11] "addr:flats" "addr:full"
## [13] "addr:hamlet" "addr:housename"
## [15] "addr:housenumber" "addr:inclusion"
## [17] "addr:interpolation" "addr:place"
## [19] "addr:postbox" "addr:postcode"
Cada feature contiene una lista de tags.
Puede acceder al vector de tags usando la función
available_tags(). Por ejemplo, puede acceder a la lista de
amenity así:
available_tags("amenity") %>% head(20)## [1] "animal_boarding" "animal_breeding" "animal_shelter"
## [4] "animal_training" "arts_centre" "atm"
## [7] "baby_hatch" "baking_oven" "bank"
## [10] "bar" "bbq" "bench"
## [13] "bicycle_parking" "bicycle_rental" "bicycle_repair_station"
## [16] "biergarten" "boat_rental" "boat_sharing"
## [19] "brothel" "bureau_de_change"
Para descargar features desde OSM, primero
debe definir un espacio geográfico y obtener la caja de de coordenadas
que lo contiene:
## obtener la caja de coordenada que contiene el polígono de Bogotá
opq(bbox = getbb("Bogotá Colombia"))## $bbox
## [1] "4.4711754,-74.2235137,4.8331695,-74.0102577"
##
## $prefix
## [1] "[out:xml][timeout:25];\n(\n"
##
## $suffix
## [1] ");\n(._;>;);\nout body;"
##
## $features
## NULL
##
## $osm_types
## [1] "node" "way" "relation"
##
## attr(,"class")
## [1] "list" "overpass_query"
## attr(,"nodes_only")
## [1] FALSE
Puede almacenar un objeto osm de clase list
y overpass_query con las estaciones de autobús para la
ciudad de Bogotá:
## objeto osm
osm = opq(bbox = getbb("Bogotá Colombia")) %>%
add_osm_feature(key="amenity" , value="bus_station")
class(osm)## [1] "list" "overpass_query"
Puede acceder a los Simple Features del objeto
osm usando la función osmdata_sf(). El objeto
contiene una lista de objetos con los puntos, líneas y polígonos
disponibles.
## extraer Simple Features Collection
osm_sf = osm %>% osmdata_sf()
osm_sf## Object of class 'osmdata' with:
## $bbox : 4.4711754,-74.2235137,4.8331695,-74.0102577
## $overpass_call : The call submitted to the overpass API
## $meta : metadata including timestamp and version numbers
## $osm_points : 'sf' Simple Features Collection with 2708 points
## $osm_lines : 'sf' Simple Features Collection with 102 linestrings
## $osm_polygons : 'sf' Simple Features Collection with 381 polygons
## $osm_multilines : NULL
## $osm_multipolygons : 'sf' Simple Features Collection with 102 multipolygons
Puede acceder a los elementos de la lista y crear un objeto de clase
sf:
## Obtener un objeto sf
bus_station = osm_sf$osm_points %>% select(osm_id,amenity)
bus_station## Simple feature collection with 2708 features and 2 fields
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: -74.22302 ymin: 4.530763 xmax: -74.03537 ymax: 4.812363
## Geodetic CRS: WGS 84
## First 10 features:
## osm_id amenity geometry
## 261630505 261630505 <NA> POINT (-74.14692 4.631477)
## 261630506 261630506 <NA> POINT (-74.14694 4.631353)
## 261630508 261630508 <NA> POINT (-74.14423 4.631102)
## 261630509 261630509 <NA> POINT (-74.14422 4.631225)
## 266537186 266537186 <NA> POINT (-74.11613 4.654154)
## 266537187 266537187 <NA> POINT (-74.11531 4.652883)
## 266537189 266537189 <NA> POINT (-74.11439 4.652935)
## 266537456 266537456 <NA> POINT (-74.11531 4.65398)
## 266537457 266537457 <NA> POINT (-74.11507 4.65568)
## 292774696 292774696 <NA> POINT (-74.14562 4.631682)
Visualiza el objeto bus_station:
## Pintar las estaciones de autobus
leaflet() %>% addTiles() %>% addCircleMarkers(data=bus_station , col="red")Aquí puede acceder a los datos del Censo Nacional de Población de 2018 para Colombia.
## censo data
browseURL("https://microdatos.dane.gov.co//catalog/643/get_microdata")
##=== variables ===##
## COD_DANE_ANM: Codigo DANE de manzana
## UA_CLASE: ID
## COD_ENCUESTAS: ID de encuesta
## U_VIVIENDA: ID de vivienda
## H_NRO_CUARTOS: Número de cuartos en total
## HA_TOT_PER: Total personas en el hogar
## V_TOT_HOG: Total de hogares en la vivienda
## VA1_ESTRATO: Estrato de la vivienda (según servicio de energía)
##=== load data ===##
## unzip file
unzip(zipfile="input/11_BOGOTA_CSV.zip" , exdir="input/." , overwrite=T)
## data manzanas
mgn <- import("input/censo_2018/CNPV2018_MGN_A2_11.CSV")
colnames(mgn)
mgn <- mgn %>% select(COD_DANE_ANM,UA_CLASE,COD_ENCUESTAS,U_VIVIENDA)
## data hogar
hog <- import("input/censo_2018/CNPV2018_2HOG_A2_11.CSV")
colnames(hog)
hog <- hog %>% select(UA_CLASE,COD_ENCUESTAS,U_VIVIENDA,H_NROHOG,H_NRO_CUARTOS,HA_TOT_PER)
## data vivienda
viv <- import("input/censo_2018/CNPV2018_1VIV_A2_11.CSV")
colnames(viv)
viv <- viv %>% select(COD_ENCUESTAS,UA_CLASE,U_VIVIENDA,V_TOT_HOG,VA1_ESTRATO)
## join hogar-vivienda
viv_hog <- left_join(hog,viv,by=c("COD_ENCUESTAS","U_VIVIENDA","UA_CLASE"))
## joing mnz-hogar-vivienda
viv_hog_mgn <- left_join(viv_hog,mgn,by=c("UA_CLASE","COD_ENCUESTAS","U_VIVIENDA"))
##=== collapse data ===##
db <- viv_hog_mgn %>%
group_by(COD_DANE_ANM) %>%
summarise(med_H_NRO_CUARTOS=median(H_NRO_CUARTOS,na.rm=T),
sum_HA_TOT_PER=sum(HA_TOT_PER,na.rm=T),
med_V_TOT_HOG=median(V_TOT_HOG,na.rm=T),
med_VA1_ESTRATO=median(VA1_ESTRATO,na.rm=T))
## export data
export(db,"output/mnz_censo_2018.rds")
## delete files
unlink("input/censo_2018/",recursive = T)
unlink("input/__MACOSX/",recursive = T)Aquí puede acceder al Marco Geoestadístico Nacional para Colombia.
## MGN data
browseURL("https://geoportal.dane.gov.co/servicios/descarga-y-metadatos/descarga-mgn-marco-geoestadistico-nacional/")
## load data
mgn <- st_read("input/mgn/MGN_URB_MANZANA.shp")## Reading layer `MGN_URB_MANZANA' from data source
## `/Users/eduard-martinez/Dropbox/teaching/bid-data-real-state/2023-01/clase-03/public/input/mgn/MGN_URB_MANZANA.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 43438 features and 13 fields
## Geometry type: POLYGON
## Dimension: XY
## Bounding box: xmin: -74.36463 ymin: 3.888847 xmax: -74.00074 ymax: 4.833006
## Geodetic CRS: WGS 84
mgn %>% head()## Simple feature collection with 6 features and 13 fields
## Geometry type: POLYGON
## Dimension: XY
## Bounding box: xmin: -74.13846 ymin: 4.661075 xmax: -74.11615 ymax: 4.683528
## Geodetic CRS: WGS 84
## DPTO_CCDGO MPIO_CCDGO CLAS_CCDGO SETR_CCDGO SECR_CCDGO CPOB_CCDGO SETU_CCDGO
## 1 11 11001 1 000 00 11001000 6402
## 2 11 11001 1 000 00 11001000 6402
## 3 11 11001 1 000 00 11001000 6403
## 4 11 11001 1 000 00 11001000 6403
## 5 11 11001 1 000 00 11001000 6302
## 6 11 11001 1 000 00 11001000 6302
## SECU_CCDGO MANZ_CCDGO MANZ_CCNCT MANZ_CAG Shape_Leng
## 1 03 23 1100110000000064020323 145420 0.006490357
## 2 09 19 1100110000000064020919 145102 0.004081371
## 3 02 15 1100110000000064030215 144714 0.002126900
## 4 02 17 1100110000000064030217 144665 0.001641266
## 5 04 10 1100110000000063020410 147399 0.002941721
## 6 02 03 1100110000000063020203 147590 0.003051028
## Shape_Area geometry
## 1 1.413612e-06 POLYGON ((-74.1378 4.677946...
## 2 9.998505e-07 POLYGON ((-74.13436 4.67957...
## 3 2.572230e-07 POLYGON ((-74.13465 4.68263...
## 4 1.500361e-07 POLYGON ((-74.1348 4.683005...
## 5 4.549105e-07 POLYGON ((-74.1215 4.663893...
## 6 6.541770e-07 POLYGON ((-74.11667 4.66107...
table(mgn$CLAS_CCDGO)##
## 1 2
## 43361 77
mgn <- mgn %>% subset(CLAS_CCDGO==1)
mgn <- mgn %>% select(MANZ_CCNCT)## load data censo 2018
censo <- import("output/mnz_censo_2018.rds")
censo## # A tibble: 38,567 × 5
## COD_DANE_ANM med_H_NRO_CUARTOS sum_HA_TOT_PER med_V_TOT_HOG med_V…¹
## <chr> <dbl> <int> <dbl> <dbl>
## 1 1100110000000000000000 NA 0 NA NA
## 2 1100110000000011010101 6 13 1 2
## 3 1100110000000011010102 3 442 1 2
## 4 1100110000000011010103 4 787 1 2
## 5 1100110000000011010104 4 2511 1 2
## 6 1100110000000011010105 2 224 1 2
## 7 1100110000000011010106 3 223 1 2
## 8 1100110000000011010107 3 279 1 2
## 9 1100110000000011010108 3 121 1 2
## 10 1100110000000011020101 4 166 1 2
## # … with 38,557 more rows, and abbreviated variable name ¹med_VA1_ESTRATO
## agregar a manzanas de mgn info de censo
mnz_censo <- left_join(mgn,censo,by=c("MANZ_CCNCT"="COD_DANE_ANM"))
mnz_censo## Simple feature collection with 43361 features and 5 fields
## Geometry type: POLYGON
## Dimension: XY
## Bounding box: xmin: -74.22353 ymin: 4.46292 xmax: -74.00074 ymax: 4.833006
## Geodetic CRS: WGS 84
## First 10 features:
## MANZ_CCNCT med_H_NRO_CUARTOS sum_HA_TOT_PER med_V_TOT_HOG
## 1 1100110000000064020323 4 465 1
## 2 1100110000000064020919 4 656 1
## 3 1100110000000064030215 NA NA NA
## 4 1100110000000064030217 4 93 1
## 5 1100110000000063020410 6 55 1
## 6 1100110000000063020203 4 538 1
## 7 1100110000000063020239 5 175 1
## 8 1100110000000063020517 4 647 1
## 9 1100110000000063020513 4 439 1
## 10 1100110000000063120117 5 104 1
## med_VA1_ESTRATO geometry
## 1 3 POLYGON ((-74.1378 4.677946...
## 2 3 POLYGON ((-74.13436 4.67957...
## 3 NA POLYGON ((-74.13465 4.68263...
## 4 3 POLYGON ((-74.1348 4.683005...
## 5 4 POLYGON ((-74.1215 4.663893...
## 6 4 POLYGON ((-74.11667 4.66107...
## 7 4 POLYGON ((-74.11826 4.66223...
## 8 4 POLYGON ((-74.11716 4.66048...
## 9 4 POLYGON ((-74.11905 4.66280...
## 10 4 POLYGON ((-74.1238 4.668334...
## visualizar estratos
ggplot() + geom_sf(data=mnz_censo , aes(fill=as.factor(med_VA1_ESTRATO)) , col=NA)## visualizar densidad poblacional
ggplot() + geom_sf(data=mnz_censo , aes(fill=sum_HA_TOT_PER) , col=NA) +
scale_fill_viridis(option="inferno")## exportar sf de manzanas con info del censo
export(mnz_censo,"output/mgn_censo_2018.rds")Lovelace, R., Nowosad, J., & Muenchow, J. (2019). Geocomputation with R. [Ver aquí]
Bivand, R. S., Pebesma, E. J., Gómez-Rubio, V., & Pebesma, E. J. (2013). Applied spatial data analysis with R. [Ver aquí]